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ABS TRACT 


The partial differential equation for heat diffusion is 
numerically integrated by the Runga-Kutta method. Solutions 
are obtained for the diurnal tepperature variation with a 
bounded coefficient of eddy diffusivity which varies with the 
lapse rate and with height. The surface wave i8 represented 
by the sum of a diurnal and a semidiurnal harmonic wave. The 
results may be interpreted to apply over a fairly broad range 
of diffusivity with height. With appropriate choices of the 
various parameters,reasonably good agreement is obtained be- 
tween theoretieal and observational values of amplitude and phase 


lag as functions of height and time. 
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Le Introduction, 

Many efforts have been made by meteorologists to solve 
the heat diffusion problem in the atmosphere, in particular the 
diurnal temperature wave has received much attention. In every 
attempt the concept of eddy diffusivity has been retained, the 
various functional forms of eddy diffusivity give different 
results. Most analyses are based upon the solution of the Taylor 
heat diffusion equation, which states that the turbulent flux of 


heat is proportional to the gradient of potential temperature, Namely: 


o fe) Pe 
mB = 2 (Kk 28 (1) 


aC rer 

Here t represents the time, Z the height, and 8G mc) the potential 
temperature deviation from a constant (mean) value. 

The simplest approach is based upon the assumption of a 
constant diffusion coefficient K. Aceording to this solution the 
amplitude of the diurnal temperature wave decreases exponentially 
and phase increases linearly with height. For a variable K 
several solutions have been presented. Sutton (1) in his compre- 
hensive survey of daily temperature variation has considered K <= Kies 
where m is the stability parameter having magnitude Ogd@m<\ . 
His final solution is in terms of Bessel functions, It may be 
inferred from his ej ka at sufficiently great heights, the 
diurnal variation becomes periodic with a pnase lag which is 
proportional to 2 » 80 that the exponent of z is now as 
low as 0.5 (for m=1). ftowever, agreement with observation 18 not as 


good as might be expected in view of the additional complexity of 


the solution. 





Best and Kohler (2) have made a similar approach to the 


problem. Haurwitz (3) has considered K as: Ke Ko + K,z , where 


1 


Ko and Ry are constants. By proper selection of Ko and K 


1, 
Haurwitz obtains good agreement with data taken at Leafield, 
England, up to 10 meters. In his results, the amplitude falls 

off more rapidly at low levels than in K Constant Theory. How- 
ever, it is greater at high levels where K is larger. The phase 
lag is larger near the ground where K is small but becomes smaller 
at higher levels where K is larger. Poppendiek (4) has obtained 

a solution, for K varying sinusoidally with time and linearly 
with height. The selution is very complicated, consisting of an 
infinite trigonometric series having coefficients determined by 
difficult integrations of functions of the real and imaginary 
parts of Hankel functions of complex argument. No numerical values 
are given, and the results are not suitable for practical purposes. 
Even for the case of K independent of time but varying as a powln 
of height, the solution is in terms of Bessel functions of order 
depending on this pewer. Staley (5) has investigated the problem 
with K increasing with height but bounded. In several respects 
the results show better agreement with observations than previous 
solutions for the coefficient unbounded. A numerical solution 
with a diffusion coefficient which varies periodically with time 
and exponentially with height has been obtained by Haltiner (7). 
In general the functional forms of the diffusion coefficient have 
not included such parameters as the roughness coefficient, gradient 


wind, lapse rate etc., which are generally believed to be closely 





related to the diffusion process. It is the purpose of this investiga- 
tion to attempt to determine a diffusion coefficient in which the time 
variation is related to the stability. Briefly, the present analysis 

consists of the selection of proper functional form for the coefficient 
of eddy diffusivity from observed cata, followed by numerical solution 


of the Taylor heat diffusion equation for this value of K. 





2. Heat Diffusivity Coefficient. 

From physical considerations, there is no reason to expect 
an ever increasing K. Hence a bounded function of K appears to 
be a logical choice. Also, many studies on the variation of K 
with various meteorological parameters give clues to the form- 
ation of a suitable function. Sutton (1) has remarked that " the 
relatively small change of amplitude of the diurnal temperature 
wave above 10 meters in sunimer compared with that in winter 
indicates that a more effective mechanism for the upward transfer 
of heat is at work in the summer than in the winter''. Also 
Haltiner (7) suggests that K should have diurnal variation. 

It is possible that both the seasonal and diurnal variation in 
diffusion may reasonably be accounted for by making K a suitable 
function of lapse rate. Moreover from the observations made at the 
micrometeorological tower at the University of Washington, Fig. (1), 
which shows K as a function of time at three heights, it is evident 

© 
that K at a given level increases rapidly by one or more orders of 
magnitude as the lapse rate changes from sub-adiabatic to super- 
adiabatic Pilar wont Wee . In the evening when the lapse rate changes 
from super-adiabatic to slightly sub-adiabatic K drops rapidly by one 
or more order of magnitude at any given level. Further from Fig.(2),_ 
which shows K as a function of lapse rate, it is evident that at lower 
levels the variation in K is less than at higher levels. On further 
investigation of the data collected at the University of Texas micro- 
meteorological tower (see Table 1) it has been found that the variation 


in K at the lower levels is an order of magnitude (factor of ten) while 








20 24 og o8 12. \G 219 
Time (Local) 


Coefficient of eddy conductivity as a function of time at 
different heights; a: 0.75 m at University of Washington micro- 
meteorological mast and tower, represents mean over 23 days of 
Aug. 1950, chosen for having most nearly sinusoidal temperature 
variations; b and c : for 15 and 35 m, respectively ( after Jehn 
and Gerherdt, 1950) , 3 - 4 Aug 1948, at Manor Texas, Dashed 
lines represent unreliability. All values obtained by methods 


based on energy continuity at earths gurface. 





Veeck oma epee oi eoean aa 


__— Te ow 





ee eT ES he Y Ay, OE Le aR GE a el OS, a, es Oo Fe er tel Mh Awe 








ST Lt OT a eto BO Rr Lees, 


26 
%e.- 





7000 = 


TO00°- 


£000 °- 


SRS ee Per OY RD LIES aE TN EE, OE ER Pe & Uwe Ar,’ 
Zz. : 
Be 


ST on OD Gaaicie ia Nae Ee 


| 
| 
| 
| 
: 
~<a 
| 
| 


PsOL ~ 8 01 X 0°8 
| 901 aed cOl X O'F 
| Ot ce ee cOT * L°2 
| got Soe! c0l X 0°2 
| co x o° gOT ¥ Z°T 
| cOT X 0°S eats G99 
i 40 
2. __ (82999H) SE 
cOl * 9°T yOl x 0°9 
401 * 7°9 cOT X 9°T 
gOT X cil cOT X S°2 
Eola xa i cOl X 9°T 
cot X 0°2 OI X 9°S 
3) a | y01 * 7°9 
7 
a >| 
(s2332H) CT 


ee ee oe ee et 





A Ay tate AK OOS Sy tS, a, pee teen, Sma atten 








TBO 
a 


Se 0 et eh Chee HE Te A Ng te Re ae a a eee een Daly tae at te ee Fie ae ge Catan, ta a LR pene, —- Pn, 0 oe 


T eTqeL 


cOl X B°2 
yO ¥ 9°2 
cOT ¥ O'T 
cOl ¥ 0'T 
cOl X 9°2 
cOT X 9°2 


TBO 
| 


y0l * 7S 


yOT X 2°S 


yOl ¥ S°P 


gol * O'S 


yOl X 0°9 


Gute 1.9 





ane Dat ree pK) Oe en - Os oe ee 


eed 


701 X S°9 €000° 
cOl ¥ S°Z £000° 
cOT XE 000° 
cOl % 2 000° 
col * £ Z000° 
oOT X 2 Z000° 
ow Bac - 
(s1923y) ¢Z 
cOl * 8 S000’ 
cOT X ZT OT 00° 
zOl x $°9 SE00° 
201 x L S100" 
mol ee eT 0200° 
70 * cl ¢200° 
a er 
(st37aW) ¢ 





Oe ow eee er me a =, 


Ke a al 
¢ 
& 
€ 


0080 


ow- A Ue . 


000 


OOVC 


000cd 


O009T 


OO0cT 


auy I, 


0080 


0070 


002 


ee ee ee i el De eh le te eel et eel See sstatl peeps an ee “Be a 8 Ca ee ee line tg ES Se fem © 


0002 


OO9T 


O0cT 


SD oo rt ee ES TYE OE 


ouTL | 


FN at AO Cee 





—N 


{x4 


CORO iret 80 es eS tS, res Ruane air a ik lle LE PELE OLDE CTL LE ETD VOD CORRS Tie IIE ty ar 


+ 


\ 


s 
\ 


» 





eS 
10 





+101 


nah An a wR HPO apnre oats nee ree Mnereeneeal 
+ » Oo i 


+ O00} 


om ¢ oo | 


-'O| 


Co 
\o 


cu) 


woe 


a, Lele 





at higher levels K varies by two orders of magnitude. Therefore 
any function which represents K should show the above character- 


istics. As a preliminary study the following functions have been 


investigated: 
B: 8 Soe 
R= oli-Gr vase) 3 |[i-ay 78" ] @) 
_a i. 28 
Ko Q, ji a, e° zy) (3) 
a6 Qa z* 2° 
K = ks ae ro 
aia, 28) oo BS (@ 
ce a} Ss. z 2 ]L- aye “f] a 


Here z represents the height, 0 # is the potential temperature 
deviation frem mean value and a. ay » a, » ay» a. are contants. 

In equation (2) the second factor ( t— a, es ) 1s obvieusly 
bounded; and the other factor (i- (a, +322 eS) changes K as lapse 
rate changes. Furthermore, the variable z is also involved in this 
factor. This is necessary, because at great heights the change in 
lapse rate becomes very small and the factor z is needed to amplify 
the effect of lapse rate at these elevations. When z = 0, oe is 
maximum, then K reduces to K a (a,-0, 28 ) ('— ag ) for zs oO » 220 & 
K=Q,In quality this equation possesses all characteristics mentioned 
above but the quantitative agreement with the observed data is not 
sufficiently good. In the second representation (3) K is bounded 
and it increases with height to a certain level which may be deter- 


mined by the proper selection of constants. There is an appreciable 
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difference in day-time and night-time, but the night-time values are lower 
at high elevations tiaalag lower levels. This function is very sensitive 
to changes of 23 and in magnitude of z. The square of the quantity 
was intended to avoid negative values of K. 
Form (4) is somewhat simpler than equation (3). In this equation also, 
K is bounded and has much the same characteristics as equation (3), but the 
variation in the value of K when lapse rate changes sign is less than the 
former. Diffusion coefficient (5) gives a reasonably good fit with the 
observed data although the night-time values are higher than observed. 
The difference between night and day-time values is correctly maintained. 
For the present analysis, equation (5) has been selected as the best 


representation for K of this group. The following numerical values are 


chosen as being representative: 


4)—2 5 x 109 cm? sec~! 5 a9 2 .2x 10°? 
a, = .99 ; @ = .2x 1073 


The choice of a, gives a hundred-fold variation of K with height. The 
value of a, controls the overall magnitude of K: 

For z=0 ,K 5x 102 cm? sec”! 

For z=00 , K 5x 1092 cm? secml 
The value of ao always keeps the value of K positive by limiting the value 
of ao 2B < for all elevations and typical lapse rates. The para- 


meter a, controls the rate of decrease of K with height. 





3. Solution of Taylor Heat Diffusion Equation . 


The associated boundary conditions of the heat diffusion equation 


(1) are: 
= o forze# © and allt (6) 
G=6(t) for z = 0 (7) 


In the analysis done by Haltiner (7) 6, (t) has been represented 
by two trigonometric terms which is normally a reasonable good approx- 
imation to the observed data. Thus the following form is assumed for: 


0 = %¢ (Sar we + ¢ Se 2 wt ) (8) 


ar -] 
Here W= 864.00 8e¢ 3 a. = .087 which for convenience keeps the 


magnitude of 8, S°!l 3; c¢ # 0.3 

In order to scale the equation for computational purposes, let 
t = 3600Cand z = 100qo (9) 
Here q is a constant, while o and € are the new variables, with t 
in seconds, the units of C 7 hours. For z in cm, o will be units 


of q-meters, ce in meters when qexl, in 2 meters unit when q = 2, 


ete., with the transformation (9) equations (1) and (4) become: 


ae _ 26 (Kae 4 BB ae (10) 
at q Poon oo AS 
— A, o~ 
crafter eyicme 7] ay 
The form of (8) remains the same except that wnowust now be taken 
as —s . The boundary conditions (6) and (7) remain identical in 


form. It will be worth mentioning here that two-meter level may be 


considered as the lowest level in order to neglect the surface effect. 
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4. Finite Difference Equations. 


In order to obtain a numerical solution of the problem the derivatives 
of §@ with respect to so in equation (10) and (11) are replaced by appro- 
priate finite difference forms. The problem is then reduced to one of 
solving a system of linear algebraic equations in the values of @ over a 


grid of points covering the desired range of time and height. We obtain: 


are bass = [KGa 28 +. y+ h os. i? ®...) (12) 


i =[1~ % oe Q,¢ “n(e,,- 8 a sI!- uae) “7 (13) 


A 


&) =|- oy (Be 4,72 82+ Be -,) . O24 Gre as _ i: * li 


[ — A a, i A (0.7 ra [os Age i (14) 
2 


flere the subscript i designates the ith level of the vertical grid at 
height ih, where h is the vertical distance between the points of the grid 
in units of q-meters. Substituting the values of (13) and (14) for i= 1, 2,--- 
O08 
in the right hand side of (12) the respective values of Bye = i= 1, 2322. 
are obtained. Now the partial differential equation (1) has now been reduced 
to a system of ordinary differential equations, which will be integrated 


numerically by the Runga - Kutta method. 


Li 





a Computation. 

For this particular problem q is selected as q = 10>, This 
means the units of o are in thousands of meters. In order to 
achieve reasonable accuracy a space interval of ten meters is 
selected in the vertical. This will be an appropriate time to 
mention that there will be less accuracy in the results below ten 
meters. All the figures below this level are approximated linearly. 
A fairly detailed picture of> the vertical structure is obtained 
assuming that the upper boundary condition (6) applies at the 
height of 290 meters. Therefore 9,4 = QO . Hence the system 
(12) to be solved consists of 28 simultaneous ordinary differ- 
ential equations, together with equation (8) for @, . A numerical 
solution of this system is obtained by the Runga - Kutta method on 
a National Cash Register 102A electronic computer. The choice of 
the appropriate time interval over which the integration is to be 
carried out depends intimately upon the size of the eddy diffusivity. 
As the latter is increased, AC must be decreased in order to 
maintain sufficient accuracy. With the value of a, po 10° » the 
time interval of the integration period was prohibitively small for 
the type of computer being used, i e., thousands of hours of machine 
time would be required to complete one full cyele. Hence a smaller 
value of a,= 1,92x10° was used, which permitted a time interval of 


1 
1/16 hours. 
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6. Results. 
5 


The reduction of ay from 5 x 10° to 1.92 x 107 wedeeas the 
coefficient of diffusivity which, in turn, reduces the amplitude 
and increases the phase leg with elevation. In Fig. (3) relative 
amplitude and phase lag in minutes are plotted against height in 
meters. It may be ebserved that the amplitude at ten meters is 
approximately one half of the surface value and the phase lag at 
this level is 7& minutes, These values are not uncommon in winter. 
These values may be compared with those observed at Porton in 
December on clear days. For instance at 17 meters at Porton 
the phase lag is 72 minutes and theory gives 63 minutes. In 
Fig. (3) for comparative purposes, the amplitude reduction and 
phase lag observed at the University of Washington micrometeor~- 
logical mast and tower is given in the inset, Notice that in 
the theoretical results amplitude falls off more rapidly at lower 
levels where K is small and the change of phase is smaller than at upper 
levels. The lower values of amplitude and larger phase lags 
compared with observed data once again indicative of the small 
value of K. The value of K in the present analysis at the surface 
is Ka 1.95 x 10% and at the topmost layer, K = 1.95 x ite. These 
are very small values particularly for clear days in August on the 
land. 

As givem in Table (1), the magnitude of lapse rate (in deg 
per cm) at lower levels, considering above the two meter level, is 
of the order, ‘hae and between ten and forty meters, this order is 


of ior, In general the lapse rate continues to decrease with elevation. 
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From the cemputation, the values at lewer levels (abeut LU meters) 
and upper levels are of order iow and to respectively which are 
about one-hundredth of the observed values. This decrease in lapse 


i 
rate may be overcome by selecting a, in equation (5) hundred times 


2 
large. In equation (12) instead of reducing a, which decreases the 


J 
magnitude of the coefficient of diffusivity, the factor q may be 
increased to achieve computational stability in the computer. The 
increase in g increases the , finite vertical distance h over which 
the devivatives are evaluated. This implies that the accuracy of 
the finite- difference approximation will decrease. 

As an example of this technique, let q be increased by a factor 
of 10, then K will be equal to 1.95 x oan which is a reasonable 
value for the diffusion coefficient. Then the finite vertical dist- 
ance h will be 100 meters. In this case Fig. (3) may be considered 
as plotted with height from zero to 3000 meters. This interpretation 
of the results may be of value if the temperature changes at higher 
elevations are of primary interest. 

Fig. (3) also shows the amplitude reduction and phase lag for 
the minimum value of the temperature at various elevations. The 
amplitude reduction for the minimum shows greater reduction than the 
maximum, which apparently corresponds to the lower nocturnal 
diffusivity. There is an appreciable difference between the amplitude 
reduction of the maximum and that of the minimum. For example, at the 
40 meters level, the maximum is about 26% of the surface value of 9 , 
while the minimum is about 19% of the surface minimum. At much lower 

and much higher levels, the difference between amplitude reduction 
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of the maximum and minimum is less. The phase lag increases with 
altitude for both maximum and minimum, but the rate of increase 
decreases with height, because of the inerease of the coefficient 
of diffusivity. The phase lag of the minimum is somewhat less 
than that of the maximum which does not seem to be consistent 
with greater amplitude reduction of minima compared to the maxima. 
Perhaps these differences would vary in succeeding maxima and 
minima a8 a greater time elapsed from the initial conditions. 

Some comparisons between the observed data and theoretical 
results may be made by a proper choice of q which indicates the 
appropriate diffusion coefficient. Table (2) shows the amplitude 
reduction and phase lag from observation taken at.Leafield. The 
agreement here is fair at elevations above 60 meters. The choice 
of qa6é gives the finite difference interval of 60 meters. Table 
(3) shows some observation taken at the University of Washington, 
Seattle. In this case the finite difference interval is 160 
meters, but observations are below this height. However, as a 
first approximation, the amplitude reduction and phase lag are 
interpolated between 0 and 160 meters and tabulated. The results 
fit well, but the accuracy of the interpolation is doubtful. 
Similarly Table (4) compares the theory with observed data at 
Eiffel Tower. Here q is very large giving rise to finite diff- 
erence of h # 600 meters, whereas the heights considered at the 
Eiffel Tower are below 300 meters. Yet the comparisen with the 
theory is fair, In Table (5) comparisen is made with the results 
taken at Porton. In this case h «= 20 meters which is not so great 


as to offset the results of computation. The theory fits well with 
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the observed data. Table (6) compares the amplitude reduction in summer 
and winter at Leafield with theory considering q = 25 and 6 respectively, 
this gives h = 250 and 60 meters respectively. The values of 9 considered 
here are proportional to the ten-meter level. The magnitude of amplitude 


reduction in theory and observation are very close to each other. 
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Table.2. 
Amplitude reduction (A.R.) and phase lag (P.L.) taken for the 


case q~ 6 versas observation taken at Leafield, England. 








Theory Leafield 
ALR. EL. Height AVR. P.L. 
iS BS 25 meters 90 - hours. 
68 42 30.=«O"" = 20! 
_ - 50 " 8h = 11 
De 3 60 " 7 1.5 «(8 
vs 1.75 90 " Z 1.66 " 
AL To HOO aie | a ee 
od 2.3 120 " = _ 1! 
Tae. 36 : 


Amplitude reduction and phase lag taken from case q - 16 versus 


observation at the University of Washineton, Seattle. 





Theory Observation 
ies al Height Aah er ee 
POL 12 min 15 meters OL 16 min 
ee) ome 510 86 if 
eo ae 45 i Roi: ao) st 
73 5S ie 60 = i 
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Table.4. 


Amplitude reduction from Eiffel tower versus theoretical values 











Tor = 60, 
Theory Observation 
Bett. Height(meters) Summer Winter 
Ah, Lo 88 83 
68 300 Plas ~60 
Table. 5. 


Amplitude reduction and phase lag taken for the case q = 2 versus 


obsertion taken at Porton, England. 











Theory Porton 
Boots Pas. Heicht A.R. Pele 
AS ced 7 meters Be is, 95 hours 
hs 05 7 65 20 
Table,6, 


Amplitude reduction from Leafield (England) versus theoretical values. 





Height in meters 10 2 50 16: 100 


Relative amplitude. 


(a) June. POU Oeod 0.84 One Oi 
(b) December. 1.00 Or7L Ope Oru. 0.40 
Theory: q = 25 1.00 0,90 0.84 0.78 0.64 

gens 1.00 0.64 0.58 0.47 0.42 
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7. Conclusions. 

A numerical solution has been presented for the diurnal temper- 
ature varlation with a coefficient of eddy diffusivity which is a 
function of height and lapse rate. By suitable selection of several 
parameters, reasonably good agreement has been obtained between 
theory and observation. 

Improvement in the results may possibly be achieved by improving 
the functional form of diffusivity coefficient. It would be partieul- 
arly desirable to make a jepetne study of the relationship between 
eddy diffusivity and such parameters as surface roughness, wind speed, 
in addition to stability. If a suitable relationship could be obtained 


for K in terms of these parameters then the solution for the diurnal 


temperature wave would have a broader application. 
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